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Abstract. We introduce a numerical scheme to approximate a 
quasi-linear hyperbolic system which models the movement of cells 
under the influence of chemotaxis. Since we expect to find solutions 
which contain vacuum parts, we propose an upwinding scheme which 
handles properly the presence of vacuum and, besides, which gives a 
good approximation of the time asymptotic states of the system. For 
this scheme we prove some basic analytical properties and study its 
stability near some of the steady states of the system. Finally, we 
present some numerical simulations which show the dependence of the 
asymptotic behavior of the solutions upon the parameters of the system. 

AMS Primary: 65M08; Secondary: 35L60, 92B05, 92C17. 

1. Introduction 

The movement of bacteria, cells or other microorganisms under the 
effect of a chemical stimulus, represented by a chemoattractant, has been 
widely studied in mathematics in the last two decades, see [181 EH ESI E6], 
and numerous models involving partial differential equations have been pro- 
posed. The basic unknowns in these chemotactic models are the density of 
individuals and the concentrations of some chemical attract ants. One of the 
most considered models is the Patlak-Keller-Segel system pT], where the 
evolution of the density of cells is described by a parabolic equation, and 
the concentration of a chemoattractant is generally given by a parabolic or 
elliptic equation, depending on the different regimes to be described and on 
authors' choices. The behavior of this system is quite well known now: in 
the one-dimensional case, the solution is always global in time [24], while 
in two and more dimensions the solutions exist globally in time or blow up 
according to the size of the initial data, see [BJ [7] and references therein. 
However, a drawback of this model is that the diffusion leads alternatively 
to a fast dissipation or an explosive behavior, and prevents us to observe 
intermediate organized structures, like aggregation patterns. 
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In this paper, we consider a quasi-linear hyperbolic system of chemo- 
taxis introduced by Gamba et al. [11] to describe the early stages of the 
vasculogenesis process, namely the formation of blood vessels networks dur- 
ing the embryonic development. The model forms a hyperbolic-parabolic 
system for the following unknowns: the density of endothelial cells 
their momentum pu{x^t) and the concentration (j){x^t) of a chemoattractant. 
In one space dimension the system reads 

Pt + (p^)x = 0, 

{pu)t + {pu^ + P{p))^^ -apu + xp(t>x , (1.1) 
(l)t = D(l)xx + ap-b(l). 

Loosely speaking, the movement of cells is directed by the gradient of the 
chemical mediator and is slowed down by the adhesion with the substratum. 
The positive constants x a measure respectively the strength of the cells 
response to the concentration of the chemical substance and the strength of 
the friction forces. Overcrowding of cells is prevented by a phenomenological 
density dependent pressure function given by the pressure law for isentropic 
gases 

P{p)=ep\ 7>1,^>0. (1.2) 

Besides, the evolution of chemoattractant is given by a linear diffusion equa- 
tion with a source term which depends on p: the chemoattractant is released 
by the cells, diffuses in the environment and is linearly degraded. The pos- 
itive parameters D^a^b are respectively its diffusion coefficient, the produc- 
tion rate and the degradation rate and the production term is proportional 
to the cell density. 

This model was introduced to mimic the results of in vitro experi- 
ments performed by Serini et al. [28j using human vascular, endothelial 
cells. Randomly seeded on the plain gel substratum, these cells migrate and 
interact together via chemotaxis signaling and, after a while, they aggregate 
to form a network of capillaries. To be more precise, different final stages 
are observed depending on the size of the initial density of cells. For low 
densities, only isolated, disconnected clusters are formed. Increasing the 
number of cells, a sharp percolative transition occurs and a capillary-like 
network appears with a characteristic chord length independent of the size 
of the initial density. Further increase of the number of cells leads to a con- 
tinuum crossover characterized by the accumulation of additional cells on 
the network chords until the structure is no longer visible. Finally, for very 
high initial densities, a continuous carpet of cells with holes, the so called 
"Swiss cheese" configuration, is observed. From a mathematical point of 
view, emerging structured patterns, such as capillary-like networks, may be 
seen as the appearance of nonconstant asymptotic solutions with vacuum, 
namely solutions composed of regions where the density is strictly positive 
and of regions where the density (of cells) vanishes. 

To reproduce the biological setting, we consider system ( |1.1| ) on a 
bounded domain [0,L] with no-flux boundary conditions, namely, for all 
time t >0, Px{0^t) =px{L,t) = for the density, pu{0,t) =pu{L^t) = for the 
momentum and (f)x{0,t) = (f)x{L,t) = for the chemical concentration. From 
the analytical side, considering the Cauchy problem on the whole space. 



A WELL-BALANCED NUMERICAL SCHEME FOR A MODEL OF CHEMOTAXIS 3 



but in all space dimensions, it is possible to prove the global existence of 
smooth solutions, if the initial data are small perturbations of a strictly 
positive constant state, see |8||9]. For the present case of the one-dimensional 
boundary value problem, when the differential part is linearized, it is possible 
to prove the global existence and the time asymptotic decay of the solutions, 
when the initial data consist of small perturbation of stable equilibrium 
constant states, see |17j. However, the analytical study in the present setting 
is still undone and is difficult in the presence of vacuum, since the hyperbolic 
part of the model degenerates as the eigenvalues coincide when the mass 
density vanishes, but see ^1 ^6] for some rigorous results about the local 
existence of solutions for related models without chemotaxis. 

Actually, preliminary numerical simulations show that, even if we start 
from strictly positive initial data, vacuum may appear in finite time. This 
occurrence is much more physically relevant with respect to the analogous 
situation in gas dynamics. If vacuum is not expected to appear really in 
gases, it is fully relevant when dealing with the density of cells, i.e.: there 
are admissible regions without cells, and in some sense it is the main goal 
of a biologically consistent model. This is a situation somewhat similar 
to what occurs when dealing with flows in rivers with shallow water type 
equations [2J. Besides, also looking at the numerical approximation, dealing 
with vacuum needs for a special care, since we have to guarantee for the 
non negativity of the solutions This can be understood at the level of 
the associated numerical flux. A numerical flux resolves the vacuum if for 
all values of the approximate solution, it is able to generate nonnegative 
solutions with a finite speed of propagation. In this paper it will be crucial 
to use schemes with this kind of property. 

Another obstacle to a serious numerical exploration is given by the 
possible lack of proper resolution of nonconstant steady states. Actually, at 
least for not too large initial masses, solutions to system (1.1) are expected 
to stabilize for large times, around some global steady states of the system 
and, when the system reaches equilibrium, the flux is expected to vanish. 
However, most of the current schemes fail to reproduce this behavior. For in- 
stance, a classical centered discretization for the source is not precise enough 
near steady states. As well known, see for instance [25lll3], this is a usual 
problem for schemes dealing with hyperbolic problems with source. This is 
why other approaches have been introduced to balance properly the fluxes 
and the sources, so giving a more accurate approximation at equilibria. In 
the case of a semilinear model, obtained from (1.1) by neglecting the drift 
term {pu^)x and taking a linear pressure, i.e. 7 = 1, Natalini and Ribot 
p5] have proposed an Asymptotically High Order method [Ij adapted to 
the case of a system with an external source term and set on a bounded 
interval. This technique increases the accuracy of the scheme for large times 
and yields a correct asymptotic stabilization near the nonconstant equilib- 
ria. In [I2l[l3||, Gosse has studied the same problem with the well-balanced 
technique [15l[14j in the framework of finite volume schemes, obtaining sim- 
ilar results. However, both techniques are limited to diagonalizable systems 
and are difficult to extend to the quasilinear model ( |1.1| ). Let us recall also 
that, in the case of the quasilinear system with linear pressure 7 = 1 and 
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without considering the presence of vacuum, Filbet and Shu [lOj have used 
the Upwinding Sources at Interface (USI) methodology [3] to obtain a well- 
balanced scheme. Their main ingredient was the hydrostatic reconstruction, 
which was originally applied to the Saint- Venant system by Audusse et al. 
[2], implying the preservation of steady states with vanishing velocity. Ac- 
tually, no boundary conditions were considered in [lOj, and so the scheme 
was not tested against nonconstant steady states. 

The main goal of this paper is quite different, and consists in introduc- 
ing a new scheme which is able to deal both with vacuum problems and with 
problems near nonconstant steady states arising from the interaction with 
the source, which is composed by the friction and the forcing chemotactic 
terms. The strategy is to adapt the ideas in [5j, to design a scheme which 
is able to balance the numerical fluxes with the source term, considering 
new interface variables and vanishing velocities. Unlike [lOj, we treat the 
damping term through the interface variables instead of using additional, 
fractional steps to integrate it in time. We prove that our scheme preserves 
the non negativity of the density and the stationary states with a vanishing 
velocity. Finally, the scheme we consider is able to treat vacuum states at 
the free boundary, which was not considered in any of the previous works. 

Before we propose and study our scheme, we make a preliminary inves- 
tigation on nonconstant stationary solutions with vacuum for system (1.1) 
and, in particular, on single bump solutions, which are positive only on one 
connected region. We restrict ourselves mostly to the case of a quadratic 
pressure 7 = 2, which is the simplest one for finding explicit expressions. Un- 
der these two restrictions, we give a complete description of the stationary 
solutions with vacuum. We also numerically show that these solutions are 
stable and that they can be found as asymptotic states of the system (1.1), 
even in the case of strictly positive initial data. Other configurations with 
several regions of positive densities are also found numerically as asymptotic 
solutions of the system, but we cannot determine for the moment which are 
the parameters leading to one configuration or another. These results show 
that model (1.1) gives a realistic description of the vessels formations, which 
can be successfully tuned against experimental data. 

This paper is organized as follows: in Section [2j we describe the sta- 
tionary solutions in the case of 7 = 2 and of one region of positive density. 
Both cases of a lateral bump and of a centered bump are considered. Then, 
in Secti on [3| we explain the construction of the scheme we use to discretized 
system ( 1.1 ) and we prove that it preserves the nonnegativity of the density, 
and the stationary states with vanishing velocities. Finally, in Section [4} we 
present some numerical results, studying the accuracy of our scheme and 
showing that the stationary solutions of Section [2] are numerically stable. 
We also explore the dependence of asymptotic solutions on the parameters 
of the system, especially on the initial mass and on the value of the adiabatic 
coefficient 7. 



2. Stationary solutions 



In this section, we analyze the existence and the structure of stationary 
solutions to system (1.1) defined on a bounded interval [0,L] with no-flux 
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boundary conditions given in one space dimension by 

pAO,-)^Px{L,-)^0, pn(0,-) = pn(L,-) = 0, ,/>,(0,-) = </>x(i^,-) = 0. (2.3) 

Remark first that considering the evolution problem with the previous 
boundary conditions (2.3) implies that the mass of the system is constant 
in time, namely that 

M= / p(x,0)dx= / p{x,t)dx, forant>0. (2.4) 

J[0,L] J[0,L] 

Therefore, the mass M will be considered in what follows as a parameter 
which characterizes stationary solutions. 



2.1. Preliminaries. The steady states of (1.1) are the solutions of the fol- 
lowing stationary system 

Mx = 0, (2.5a) 
{pu^ + P{p))^ = -apu^XP^x, (2.5b) 
-D(j)xx^cip-H- (2.5c) 
Using the boundary condition we have that pu — Q^ and equation ( |2.5b| ) 
becomes 

P{p)x = Xp(t>x, (2.6) 
which is automatically satisfied if p and (j) are constant. More generally, for 
a pressure P{p) —ep^ with 7 > 1, we have two possible solutions for equation 
(2.6), which are 

p = 

or 

p^-\x) = ^fc^(/)(x) + K, 

where K is an integration constant. Remark that functions p and are also 
constrained to respect relation (2.5c). 

As a consequence, nonconstant equilibria can be composed of intervals 
where the density is strictly positive, which we call bumps, and of intervals 
where the density vanishes. However, we are not able to determine a priori 
the number of such intervals as a function of the data. This is why we fix 
the number of intervals with positive density and we denote it by pEN. In 
addition, we denote by xi the boundaries of these intervals, setting = 
and X2p-i — L^ so that the general form of nonconstant steady states is 

V 

k=l 

where XA is the indicator function of the interval A. On the intervals of the 
form [x2k-2^X2k-i]^ the density satisfies the relation 



1 



where (^k is the solution to 

apk-b(f)k for xE [x2/c-2,^2/c-i], A:=l,...,p, 



-h(j)k for xE [x2/c-i,X2/c], A: = l,...,p-1. 
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Finding from the above equations for a general adiabatic exponent 7 > 1 
is equivalent to solving a second order, nonlinear differential equation of the 
form 

1 

where Ai,A2,As are constants. This is a non trivial task apart from the 
case 7 = 2. 

Moreover, a unique explicit solution can be obtained only in the case 
of one interval where the density is strictly positive. Otherwise there are 
more constants to determine than available equations. More precisely, we 
need to find 4p — 2 constants coming from the integration, 2p — 2 interface 
points X and p constants K. It gives 7p — 4 parameters to determine. On 
the other hand, we have only Qp — 3 equations coming from the boundary 
conditions, the continuity of and its first and second space derivatives and 
from the total mass conservation condition. These two quantities match 
only in the case of For p>l, we would have an infinite number of 

stationary conditions depending on p—1 parameters. 

Here, we restrict ourselves to the case of only one interval where the 
density is positive. In this particular case, we can determine exactly the 
stationary solution as a function of the initial mass of the density M, the 
length of the domain L and the other parameters of the system a, ^5 7? 
a, b and D. 

2.2. Quadratic pressure 7 = 2: the case of a lateral bump. We give a 
precise characterization of the equilibria for P{p) —ep^ when we restrict our 
attention to solutions with only one bump, i.e.: p—1. One of the possible 
forms is a lateral bump, that is: there exists xG (0,L) such that the density 
satisfies 

X 




+ for xG [0,x], 
for X E (x,L], 



(2.7) 



and (j) satisfies equation (2.5c). In the following proposition, we completely 



describe these steady states, as a function of the mass M of the density 



defined in equation (2.4) 



Proposition 2.1. Consider system (2.5) on the interval [0,L] with bound- 
ary conditions (2.3) and with a quadratic pressure P{p)—ep^. // t = 

— ( — — 6)>0 and L>^^, there exists a unique, positive solution with 
D \2e / yr 

a density of the form (2.7) and of mass M . This stationary solution is 
defined by 

/ N ; ^:^(^)+^' /orxG[0,x], 
p(x) = <( 2^ (2.8a) 

0, for xe{x^L], 



and 



' 2ebK cos(Vtx) aK 
rxD cos{^/rx) rD ' 

2eK cosh(J-^(x-L)) 



X 



cosh(^ 



for X G [0,x], 
for X G (x,L]. 



(2.8b) 
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The free boundary point x is given by the only value x G ^^(7r/2,7r) 
such that 



-^ta,n(y/rx) =tanh(y -^(x — L)) 
tD V D 



and the constant K is equal to 



K- 



D 



b tan(y^x) — ^/TX 



(2.8c) 



(2.8d) 



TT 



If r<0, or r>0 but L<—=, then there is no one bump solution to the 
problem. 



Proof To p rove t he p articular form of the equihbrium (2.8), let us first 

insert (2.8a) into (2.5c). Assuming — 6)>0 and using boundary 

conditions (j)x(0^ •) ■ 



D\2e ^ 

.(L,-) = 0, the concentration (f) can be written as 
aK 

Ci cos{^/Tx) — , for X G [0, x] , 

tD 

C2COsh(w -^(x — L)), for xG (x,L]. 



Using the continuity of functions 0, (j)x and (j)xx at the interface point x, we 
obtain the values of the constants Ci and C2 , namely 

Ci = ;-cos(VTx) , 6/2 = cosh(W — (x-Lj) . 

rxD X y D 

which yield (2.8a) and ( 2.8b[ ). The remaining equation gives the relation 
(2.8c) for the location of the interface x. The parameter K can be calculated 
from the total mass density 

M= [ p{x)dx= I 
Jo Jo 



X 

Ye' 



(x)dx + Kx 



giving (2.8d). 



Now we have to determine whether there exists a solution to (2.8c) 



or not. The function /(x) 



tD 



t^in.{^/rx) — tanh I 



(x — L) is not 



defined on 



and has a positive derivative elsewhere. Its value 

/(O) is positive and the sign of f{-^ (7r + /ctt)), G Z, is the same as the sign 

of L — -j= (tt + kjv) . So, if there exists /c G N such that L > , / has exactly 

one zero in each interval of the form ((7r/2,7r) +7r{0, 1, /c — 1}). All 
these zeros are candidates to be a boundary x. However, only the smallest 
of these points, namely the one belonging to the interval ^^(7r/2,7r), guar- 
antees the positivity of the function p defined by (2.8a). Remark that in the 



TT 



case when L< there is no stationary solution of the form (2.8) 
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Finally, in the case r < 0, the same computation leads to write the 
concentration (j) under the form 



aK 



Ci cosh( v^x) — , for X G [0, x] , 

tD 

C2COsh(y — (x — L)), for x G (x,L], 



(2.9) 



and the continuity of 
interface point x 



^5 V^X5 V^XX 



at X, gives the following relation for the 



- tanh(\/— rx) =tanh \ y ^{x — L)\. 
Dr \ s D I 



(2.10) 



However, this equation has no solution for x G (0,L), since the left-hand side 
is negative and the right-hand side positive. □ 



2.3. Quadratic pressure 7 = 2 : case of a centered bump. Let us 

consider the case of a centered bump, namely the case of one interval where 
the density is strictly positive in the interior of the domain. Denoting by 
x,yG(0,L) the boundaries of this interval, with x<y, the density of the 
nonconstant steady state can be written as 



p{x) 




for X G [0,x), 
x) + for X G [x,y], 
for X G {y^L]. 



(2.11) 



In the following proposition, we describe these stationary solutions as a 
function of the mass M and of the length of the domain L. 



Proposition 2.2. Let us consider the system (2.5) set on the interval [0,L] 
with boundary conditions (2.3) and with a quadratic pressure P{p)—ep^. If 



1 



7^ 



h] > and L > —= then there exists a unique, positive solution 

2e / VT 



with a density of mass M and of the form (2.11). This solution is given by 
the following expressions: 



0, /orxG[0,x), 
p{x)={ — 0(x)+K, /orxG[x,y], 
0, forxe{y,L], 



(2.12a) 
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and 



2eK 



cosh('^ 



X 



cosh('^ 



2ehK cos(v^(x - f )) aK 
2eK ^^^Hy J^(x-^)) 



X 



for X E [0,x), 



, forx^[x,y\, 



forxe{y,L]. 



(2.12b) 



cosh(i 



The free boundary point y is such that y — L — x and so the solu- 
tion is symmetric and the boundary x is given by the only value x G 

L 7T L 7T 



2 V^'2 2^ 

tanh( 



such that 
~b . 



The constant K is equal to 



tD 

Mr 



tan(VT(x-L/2)). 



(2x-L)— -2W-^tanh(^ 



(2.12c) 



(2.12d) 



For other values of the parameters there are no solutions of this form. 

Proof. We follow the same computations as in the case of a lateral bump of 
Subsection 2^ Under the assumption r > and using boundary conditions 

for X G [0,x), 



(2.3), the concentration (f) can be written as 



Cicosh(W^ 



x) 



aK 



C2 cos(\/tx) + C3 siii{y/rx) — , for x G [x, , 

tU 



C4Cosh(W^(x-L)), 



for X G (y^L]. 



Solving the system of three equations given by the continuity of 0, (^x and (^xx 
at the interface point x (resp. y) gives an expression for the three constants 
Ci (resp. C4), C2 and C3 as a function of x (resp. y). Therefore, we obtain 
two expressions for C2 and C3, which give us two non linear equations for x 
and namely: 



sin(\/T3:^)tanh(y — x) — sin(\/T^) tanh(y —{y — L)) 



D 



v 



-^{cos{y^y) 
-cos(\/tx)). 



cos(\/T3:^)tanh(y —x) — cos(\/TI/)tanh(y —{y — L)) — 



-^(sin(VTx) 
-sin(VTy)). 
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This system can be easily rewritten as 




D 

L))sin(vr(^ 



L))cos(vr(?/-x))- 



-sin(Vr(?/-x)), 




cos(vr(?/-x)). 



(2.13) 




From ( |2.13| ) we get the relation 

h cos(v^(y-x))-l 



tanh( 



which implies that y - 



-tanh( 




Dt sin(v^(?/ — x)) 
L — x with xg(0,L/2). Inserting this relation in 



equation (2.13), we find that x satisfies equation (2.12c) 



Following the same analysis as in the proof of Prop. |2.H equation 

27r ''2i7T 
(2.12c) has a solution iff L>—= and in the case when L>—=, the only 



solution X leading to a positive density is such that ( — — 

The parameter K is computed thanks to the value of the mass M of 
the density and is given by equation ( 2.12d[ ). □ 



Remark 2.3. Computing solutions {x,y) to systems (2.13) is, in general, a 
non trivial task. In the case of one lateral bump, we found the conditions that 
guarantee the existence of x giving nonnegative density everywhere. Never- 



theless, due to the nonlinearity of equation (2.8c); we have to solve it nu 



merically to find the value of x. In the case of one centered bump, we find 
out that the solution is symmetric, which simplifies the computations. The 
conditions for the existence of solutions can be found and the explicit solu- 
tion can be calculated. However, the technique we use in that proof cannot 
be generalized to a higher number of bumps. 

Remark 2.4. Notice that the solution for the lateral bump calculated with 
L/2 and M/2 and symmetrized to be defined on the whole interval [0,L] 
is equal to the one computed for the centered bump. Our computation was 
essentially aimed to exclude other cases. 

3. Numerical approximation 
Let us now explain how to construct a reliable scheme in order to 



perform numerical simulations of system (1.1). As a standard guess, we can 



expect that solutions stabilize on steady states. This scheme will also enable 
us to test the stability of the stationary solutions computed in the previous 
section. 



System (1.1) couples equations of different natures, i.e. a quasi-linear 



system of conservation laws with sources, coupled with a linear parabolic 
equation for the evolution of the chemoattractant. The parabolic part can 
be approximated using, for example, the classical explicit-implicit Crank- 
Nicholson method. 

Now, denoting hjU= (p^puY the vector of the two unknowns, density 



and momentum, the hyperbolic part of system (1.1) can be written in the 
following form 

Ut + FiU):c = SiU), (3.14a) 
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where F is the flux function and S the source term, i.e. 

F(f)=(j;) = (^„4'*p,^)), S(U)=(l^^^J. (3.14b) 



In this section we present a finite volume scheme for (3.14) defined on a 
bounded domain [0,L] with no- flux boundary conditions (2.3). The scheme 
needs to preserve the non negativity of density and all the steady states of 
the system. 

3.1. Well-balanced scheme. According to the framework of finite volume 
schemes, we divide the interval [0,L] into N cells = [x^_i/25^i+i/2)5 cen- 
tered at nodes Xi,l<i<N. In the following, we will assume, for simplicity, 
that all the cells have the same length Ax = x^+i/2 — ^i-i/2- We consider as 



a semi-discrete approximation of the solution U of system (3.14) on cell Ci 
an approximation of the cell average of the solution at time t > 0, that is to 

say 

1 r^^+l/2 

U^{t) = -^ / U{x,t)dx. 



A general semi-discrete, finite volume scheme for ( |3.14 ) can be defined as 



AxjUi{t)+J'i+y2{t)-J'i-i/2{t) = Si{t), (3.15) 

where ^i+i/2(^) is an approximation of the flux F(U{xj^^i/2^t)) at the inter- 
face point Xi^i/2 at time t and Si{t) is an approximation of the source term 
S{U) on the cell Ci at time t. 

A classical choice is to take J^i^i/2{t) — J^{Ui{t),Ui^i{t)), where T is 
any consistent numerical flux function for the homogeneous problem 
Ut-\-F{U)x = 0- The numerical source is given as Si{t) = S{Ui{t)), where 
we discretized the derivative (j)x with a space centered formula (l)x\Ci ^ 

where (j)i±i is an approximation of function (j) at points Xi±i. 

However, it is known that this kind of approximation produces large errors 
near nonconstant steady states. 

Balancing the flux term and the source term increases significantly 
the accuracy near steady states, by imposing an exact discretization of the 
stationary solutions of the system. A possible approach is to calculate the 
flux terms J^i±i/2 in ( 3.15[ ) as a function of new interface variables U^-^i^^ i.e. 



^i+i/2 = ^(^^+1/2' ^2+1/2)' These interface variables will be made precise 
later on and their computation will take into account the balance between 
the flux term and the source. 

The technique is also to upwind the source term, defined as Si — 
^z+i/2 + ^i^i/2' ^Pi^i^ USI method [HI [271 EOj [5] . We consider 

the following ansatz: 

^.+ 1/2= T^(r.- '^^1/2= r.-^ \ ' (^'l^) 



This ansatz is motivated by an exact discretization of the stationary part 



of system (3.14a), F(U)x = S{U), using that, in the case of a stationary 



solution, the momentum pu vanishes thanks to boundary conditions. 
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Therefore, the final scheme can be written as : 



i+1/2' ^2+1/2) 



1/2! 



-1/2- 



(3.17) 



We will precise in the following subsection how to reconstruct the interface 
variables ^^^^;l/2' ^^^i^]^ contain information about the sources. This will 
be done according to the local equilibrium in order to make the scheme 



consistent with (3.14), preserving the non negativity of the density and 



preserving the steady states of (3.14) with a vanishing velocity. 



3.2. Reconstruction. In order to complete the construction of the scheme, 
we define the interface variables U^-^i^. To satisfy the well-balanced prop- 
erty and increase the accuracy of the approximation near nonconstant steady 
states, the reconstruction is obtained from the stationary system 



The system (3.18) can be rewritten in terms of the internal energy 



function e{p) which, for a pressure law of isentropic gas dynamics (1.2), is 



defined by e{p) -- 



p{p) 



*(p):=e(p) + 



We consider the function 
P{p) £7 



P 



7-1 



^ ^ for 7 > 1, 



and we divide the second equation of (3.18) by p, which leads to : 

= 0, 




+^{p)-x(i) 



-au. 



(3.19) 



We integrate now the previous system on (resp.[x^^i/25^i+i]) 
find the interface variables 



U: 



-1/2" 



^i+1/2 



resp. 



•1/2- 



^i+1/2 



and we obtain the two following equations for p^ 



i+l/2^i+l/2 



and 



g7 
7-1 



^+1/2 
,7+1 



U. 



-1/2" 



and 



i+l/2- 

(3.20a) 



i+l/2^ 



(3.20b) 



using ( |3.20a[ ), with 

£7 



£7-1 1 2 f^i+i/i 

'^iPi +xi<t^i-<l>i+i/2)-^Ui+a u{x)dx. 



For an integer 7> 1, equation ( |3.20b[ ) is a polynomial of order larger 
than two. The main difficulty in the reconstruction lies in finding its roots 
and checking that the form of U^-^i^ leads to a consistent scheme that pre- 
serves the non negativity of the density. Audusse et al. in introduced 
the hydrostatic reconstruction for shallow water equations, assuming that 
the velocity is zero at the steady states. This hypothesis simplifies equation 
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(3.20b) such that an exphcit solution can be found. This method generates 



a scheme, which is weh-balanced at equihbria with vanishing velocity. Re- 



mark that, for a system of type (1.1) set on a bounded domain with no-flux 



boundary conditions (2.3), such equilibria are the only possible stationary 
solutions of the system. 



However, in the quasilinear model of chemotaxis (1.1), the source in 



the momentum balance equation contains a damping term together with a 
chemotaxis term. The assumption i/, = at a steady state cancels the friction 
term in the reconstruction and, to satisfy the consistency property, this term 



has to be added separately in the discretization of (|3.14|), namely 





d 



-1/2 



-apiUi 



This approach was used by Filbet and Shu in [lOj in the case of a linear 
pressure function, i.e. 7 = 1. In order to include the friction term into 
the reconstruction, we define the new interface variables taking into account 
stationary solutions with constant, instead of vanishing, velocity. This allows 
to deal with the damping term together with the chemotaxis term and the 
flux term, but also to simplify equation (3.20b). Therefore, considering a 
velocity constant in space and integrating equations (3.19), we obtain the 
two following relations for uT^-^^^ and P^_^i/2- 



u. 



i+1/2 " 



and 



-1/2 



—a I 

J X 



u{x)dx + x{4^i+i/2-4^i 



(3.21a) 



(3.21b) 



Remark that for 7 > 1 the function p^^{p) is strictly increasing and contin- 
uous on [0,+oo) with a finite value at 0. So, there exists an inverse function 
which enables us to find a solution to this last equation. 



It remains now to explain how to discretize the integral in (3.21b) 
and how to find the approximation 0^+1/2- The integral in (3.21b) can be 
discretized by any consistent method, for example 

hl/2 

u{x)dx'^ —a/S.x{ui)j^^ 



-a 



I 



where X+ = max(0,X), X_ = min(0,X). The computation of ^^+1/2 is not 
completely plain. The values 0^, which approximate the function (j) at points 
Xi, are easily computed thanks to the parabolic equation for (j) and we use 
these values to calculate 0^+1/2- In order to preserve the non negativity of 
the density p, we take 0^^1/2 = niin((/)^,0^+i). Other choices, as for instance 
taking the average between and do not guarantee this property. 

In conclusion, the reconstruction of the densities becomes 



Pi^i/2^^ M [^^(pi)-aAx(i/i)+ + x(min(0i,0i+i)-0i)]+ y, 
Pt+i/2 = ( [*(Pi+i) + «Ax('^i+i)_ +x(min(0^,(/)i+i) - (/)i+i)] + ) , 

(3.22) 

where the positivity-preserving truncations guarantee the non negativity of 
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3.3. Properties of the semi-discrete scheme. In the following theorem, 
we prove some properties of the semi-discrete scheme defined by equations 
( [3l7| )-( [3l6l ) with the reconstruction ( [3^2Ta| )-( [3^ . 

Theorem 3.1. Let us consider the system ( |3.14 ) set on the interval [0,L] 
with boundary conditions ( |2.3[ ). Let T = (T^ ^T^^Y be a consistent, nu- 
merical flux preserving the non negativity of p for the homogeneous part of 
system (3.14)^ Ut^F(U)x — ^' The finite volume scheme (3.17)-(3.16) with 
the reconstruction ( |3.21a[ )-( [3.22 ).- 

(i) is consistent with (3.14) away from the vacuum, 
(a) preserves the non negativity of p, 

(Hi) preserves the steady states given by (3.19) with a vanishing velocity. 

Proof, (i) To prove the consistency of the numerical scheme with system 
( 3.14[ ), we first need to show the consistency of the flux term, i.e. 



:1,2,...,7V 



lim 



-F{U). 



Ui,Ui^i^U,Ax^O 

This is straightforward using Taylor expansions, since 

Ur^y^ = U^ + 0{Ax) and U^^y^ = U^+l + 0{Ax), 
and, therefore, 

T(u-^^/„Ul^/,)^T{Ui,Ui+i) + 0{Ax). 

The consistency finally comes from the consistency of the numerical flux T 
with the analytical flux F, i.e. ^{U.U) = F{U). 

We have now to prove the consistency of the discretization of the 
source term. To do so, we use the definition given by Perthame and Simeoni 
in [^J and we show that 



Vi = l,2,...,7V 



lim — — 

Ui,Ui^i^U, Ax^O Ax 



'^i+l/2 + '^i+l/2 



-S{U). 



Using equations (3.16|), we find that 



z+l/2^ 2+1/2" 



Pip-+l/2)-PiPi) + PiPi+l)-PiPtl/2) 



We use now Taylor expansions of ^(^^1/2)- consider the density away 
from vacuum, so that, for Ax small enough, the positivity-preserving trun- 
cations in the reconstruction ( |3.22| ) can be omitted and we obtain 

P(p7+1/2) - Pipt+l/2)-P(Pi) + Pi(-^i^i)+ + f - \<t^x,i\))Ax 



■Pipi+i)-pi+iia{ui+i)- 



X, 



^x,i+i + |</'x,i+i|))Ax + 0(Ax2 



which proves (i). 



(ii) We write the first component of scheme (3.17) : 
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To prove the conservation of the positivity of the density for our scheme, let 
us show that, whenever pi(t) vanishes, the inequality 



holds true. Considering separately the two cases 
can prove that, if Pi = 0, then 

= [^^(pi) - a(iXi)+Ax + x(min((^ 



1/2^ 

i>i > and 
H,(t>i+l)-(t>i)]- 



:0 



H+1 



we 



which implies that p[ 



-0. In the same way, we can prove that if Pi = 0, 



then PjLi/2~^- Since the numerical flux T preserves the non negativity of 
p for the homogeneous part of system ( 3.14| ), we have 
J'''iUi,Ui+i)-J^P{Ui-i,Ui)<0, 

whenever pi{t) vanishes, which completes the proof of (ii). 

(Hi). We consider a discrete version of the stationary solutions defined 
by ( |3.19[ ) with a vanishing velocity, satisfying therefore 

(3.23) 



From equations (3.22) and (3.23), we can see easily that, in that case, 
^z+i/2 ~ ^i+i/2' Stationary solutions (3.23) are preserved by scheme (3.17) 
iff 



^(^i+l/2'^i+l/2) 



^(^i-1/2 ' ^i~-l/2) ~ ^i+1/2 ~ "^t 



1/2 



i+1/2 



-1/2- 



= 0, 



which is clearly true using the consistency of T and equations (3.14b) and 
( [3l6l ). □ 

3.4. Properties of the fully discrete scheme. Let us now consider a 
time discretization of system (3.14) with a given time step At and dis- 
cretization times t^ = nAt,nGN. Using a standard approximation of the 
time derivative in scheme ( |3.17 ), the fully discrete scheme can be written as 

-1/2) ^\^i-l/2^^i-l/2)) 

At 



n+l . 



Ax 



+ 



Ax 



-1/2 



/2 



(3.24) 



where Uf- is an approximation of the solution U of system (3.14) on cell Ci 
at time t", U^J^i2 ^-re the values of the interface variables at time t", and 



5"; 







1/2" 



i+1/2 



1/2" 




P 



Pi-l/2, 



(3.25) 



We can easily see that the time integration preserves two of the properties 
proved in Theorem |3.1[ namely the consistency and the conservation of the 
stationary solutions. However, we have to find a suitable stability condition 
for the scheme, that is to say the relation between the space step and the 
time step to preserve the non negativity of the density. To establish it, we 
use the notion of invariant domain by interface given in [4J and we follow 
the proof presented in [2j. Indeed, proving directly that the scheme (3.24) 
preserves a convex domain, such as the positive half-plane, is a hard task 
since the stencil of scheme (3.24) is composed of three points. To simplify the 
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computations, we consider the weaker notion of preservation of the domain 
by interface which consists in proving two inequahties involving two points 
each. However, it is proved in [4J that these two notions are equivalent under 
a slightly more restrictive stability condition linking the time step, the space 
step and a numerical velocity. 

We first give the definition of a solver preserving the non negativity by 
interface for a fully discrete scheme, in the case of a homogeneous system. 

Definition 3.2. A solver J' = {J'^ ^J'^'^Y for the homogeneous system Ut + 
F(U)x = preserves the non negativity of p by interface with a numerical 
speed cr(Ul^^Ul^_^i) >0 if whenever the CFL stability condition 



a{Ur,Ur+i)At<Ax 



holds, we have 



p?-^{^'iur,ur+i)-p7uf)>o, 

Let us assume that the numerical flux we use for the flux discretization 
preserves non negativity by interface. The following proposition gives the 
stability condition to conserve this property in the case of system (3.14) 
with source term. 

Proposition 3.3. Let us assume that the homogeneous flux T preserves the 
non negativity of p by interface. Then the fully discrete scheme (3.24) with 
the reconstruction at interfaces given by (3.21a)-(3.22) preserves the non 
negativity of p by interface, which means than whenever the CFL stability 
condition 

^iKfi/2^K-;-t/2)^t<Ax (3.26) 

holds, we have 

At 



Proof Assuming the CFL stability condition ( |3.26 ) holds, we have 

^^+V2~ Ax V^^+l/2'^i+l/2~^ V^i+l/2'^i+l/2>'y 



(3.27) 



thanks to the Definition |3.2| of the preservation of the non negativity for the 
homogeneous system. It is equivalent to 



(3.28) 



A WELL-BALANCED NUMERICAL SCHEME FOR A MODEL OF CHEMOTAXIS 17 



From the reconstruction (3.22) we know that P^_^i/2 — — Pi+i- 
So, as long as 

At ^\ _ At 



conditions (3.28) imply equations ( |3.27 ) 



Since ^ ^(^J^i/2'^i+i/2) i = 1,2, A^, under the CFL con- 

J I Ax Ax 

dition (3.26), we have > '^i+i ^ which ends the proof. □ 



4. Numerical results 
In this section, we analyze numerically the asymptotic behavior of 



system (1.1) using the scheme introduced in the previous section. First, 
we compare three different Riemann solvers. Roe, HLL and Suliciu and we 
show that Suliciu solver is the most adapted to treat vacuum. Then, in order 
to study the accuracy of our scheme, we compare it with a standard finite 
difference method with centered in space discretization of the source term. 
We will see that our scheme captures better the interface with vacuum, 
shows less diffusion than the standard finite difference method and gives a 
proper resolution of nonconstant steady states. Therefore, for the following 
numerical simulations, we use the Suliciu relaxation solver, described in [4J 
for the system of isentropic gas dynamics, with an upwinding of the source 
term. 

This scheme is accurate enough to study numerically the stability of 
the lateral bump and the dependence of the asymptotic solutions of system 



(1.1) on some of the parameters of the system. In particular, for 7 = 2, we 
are interested in the asymptotic number of bumps for different lengths of the 
domain L and values of the chemotactic sensitivities x- We also compare 
the behavior of the system for different values of the adiabatic exponent 7 
and, finally, we study the influence of the initial mass on the structure of 
asymptotic equilibria, comparing the results for 7 = 2 and 7 = 3. 

4.1. Comparison of different solvers for the homogeneous part. In 

order to obtain a numerical approximation using a finite volume scheme we 
have to calculate numerical fluxes between control cells. This fluxes eval- 
uator is based on solving the Riemann problem at each facet in exact or 



approximate form. Hyperbolic, homogeneous part of the model (1.1) co- 
incides with the isentropic gas dynamics system for which many Riemann 
solvers are available. However, due to the presence of the source term some 
of analytical properties of solutions are modified leading to numerical diffi- 
culties that cannot be handle by all known Riemann solvers. In particular, 
occurrence of vacuum and asymptotic, nonconstant states may cause in- 
stabilities and negative values of the density. This is the reason why we 
compare different approximate Riemann solvers. Roe's method, HLL and 
Suliciu solvers and explain our choice of Suliciu as the most adapted for 
numerical analysis in the following sections. 



In the first simulation we consider system (1.1) with e — D 



L = l, x = 50, a quadratic pressure 7 = 2 and the initial data defined as 
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follows 

Po{x) = l + sin(47r|x — L/4|), (/>o = 0, uq^O. 
Figure [l] presents the (on the left) and (on the right) numerical errors 
for the three previously mentioned solvers with well-balancing of the source 
term. Behavior at two different grid sizes Ax = 0.05 and Ax = 0.01 is studied. 
The reference solution is obtained using Suliciu solver and well-balancing 
reconstruction on the fine grid with mesh size Ax = 10~^. The time step 
satisfies At = 0.9Ax/A. We observe similar behavior of the errors for all the 
solvers. They oscillate at the beginning of the evolution and stabilize with 
time. There is no significant difference between them, although the errors 
for Roe's method seem a little bigger than in the case of HLL and Suliciu. 




Figure 1. Time evolution of (on the left) and (on the 
right) numerical errors of the density p for system ( |1.1[ ) in the 
case P{p)^ep^^ with 5 = L) = a = 6 = L = 1, x = 50, approx- 
imated using a finite volume, well-balanced scheme (3.24). 
Three different Riemann solvers: Roe method (in solid blue 
line), HLL solver (in dashed green line), Suliciu relaxation 
solver (in dotted red line) are compared for Ax = 0.05 and 
Ax = 0.01 . 



In the second test, we increase the adiabatic coefficient 7 in the 
pressure function taking 7 = 3. Figure [2] presents the density profiles at 
asymptotic states containing vacuum and steep gradients near the interfaces 
between regions where the density is strictly positive and regions where 
the density vanishes. We observe that Roe's method, which is based on a 
linearization of the system, fails and produces negative values of the density 
near vacuum. We mention here that the classical Stager- Warming flux 
splitting also oscillates at the interface with vacuum. The two other solvers, 
HLL and Suliciu, are stable, they preserve non negativity of the density and 
approximate the interface with high resolution. Suliciu's approach is to use 
a relaxation scheme in which the mass conservation equation is not relaxed. 
This makes it less diffusive than the HLL solver and allows to capture 
better contact discontinuities. Moreover, in general, the HLL solver needs 
very careful wave speed estimates. However, in the case of the isentropic 
gas equations, it is not clear a priori that there is a difference between these 
two solvers and, in our tests, we don't observe any significant distinction 
between these two methods. As our model contains source terms leading to 
steep gradients of the density and appearance of regions where the density 
vanishes, we decide to use Suliciu solver, which is adapted to treat vacuum. 
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Roe 



HLL 



Suliciu 



Figure 2 
in the case P{p) 
10, M 



Density profiles for system (1.1) at time r = 5 



ep-, with £ = 1,^ = 0.1, a = 20, b 
1 approximated using a finite volume, well-balanced 



scheme (3.24). Three different approximate Riemann solvers: 
Roe method (on the left), HLL solver (in the middle), Suliciu 
relaxation solver (on the right) are compared. 



Remark 4.1. Numerical simulations indicate that none of the three approx- 
imate Riemann solvers that we studied is stable for 7 > 2 with large initial 



masses. Therefore, our numerical analysis of the system (1.1) is performed 



for initial masses small enough to assure the stability of the scheme. 

4.2. Accuracy of the numerical approximation for the source term. 



Now, we analyze how the finite volume numerical scheme ( 3.24 )-( 3.25) with 



the reconstruction (3.21a)-(3.22) manages to capture a particular asymptotic 
behavior of system (1.1). More precisely, we show that the finite volume ap- 
proach with an upwinding of the source term behaves better near noncon- 
stant steady states than a classical finite difference centered discretization. 
In this subsection, we consider the case of a lateral bump for 7 = 2. 

Figure [3] displays asymptotic density profiles at time T— 100 obtained 
with various schemes - finite volume with upwinding, finite volume without 
upwinding and finite difference without upwinding - and the exact solution 



given by (2.8). On the left, we show the computations with a space step Ax 
equal to 0.05 and on the right, equal to 0.01. The initial data are the same 
as in the previous subsection and the time step satisfies the same stability 
condition At = 0.9Ax/A. We see that the finite volume scheme with the well- 
balanced property gives clearly the most accurate location of the interface. 
Indeed, studying the different approximations of the source term, we notice 
that the centered discretization produces a bigger error at the steady state. 
Then, comparing finite difference and finite volume, we observe that the 
finite difference method is characterized by a very high numerical diffusion, 
in contrast to the finite volume approach. However, decreasing the space 
step, all the schemes converge to the reference solution, even if the finite 
volume well-balanced scheme is still the most accurate to give the correct 
location of the free boundary. 

Another important point in the numerical approximation of the solu- 



tions of system (1.1) defined on a bounded domain with no-flux boundary 
conditions is the accuracy of the velocity. In this case, the momentum 
should vanish at steady states. Figure [4] presents the asymptotic states of 
the density and the momentum obtained by two different methods : on top, 
the finite volume well-balanced scheme and on bottom, the finite difference 
scheme with a centered discretization of the source term. On the right, we 
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"0 0.1 0.2 0.3 0.4 0.5 



Ax = 0.05 



Ax = 0.01 



Figure 3. Density profiles for system (1.1) with P(p)=p^, to- 
tal mass M = l and L = l. The other parameters of the system 
are D = 0.1, a = 20, 6 = 10, x = 10. Comparison between different 
methods of approximation of a source term: (green) VW - finite 
volume method with well-balanced reconstruction; (pink) VC - fi- 
nite volume method with centered in space discretization of the 
source term, (blue) DC - finite difference scheme with centered 
in space approximation of the source. Reference solution (red) is 



given by the exact solution (2.8) 



can see the corresponding residues of the momentum. First, we observe that 
at time T= 100 the residues are less than 10~^ and decreasing, which means 
that the steady state is reached. Then, we see that the approximation of 
the density is comparable for both schemes, which was also observed in the 
previous test. However, the finite difference approach produces an error in 
the momentum profile. Its norm is close to one instead of being equal 
to zero. 




0.2 0.4 




Figure 4. On the left: Asymptotic density (blue) and momen- 
tum (green) profiles for system (1.1) in the case P{p)=£p^ and 
e = D = a = h = L = l^ X = 50 obtained by: (on the top) a finite dif- 
ference scheme with centered in space approximation of the source 
; (on the bottom) a finite volume, well-balanced scheme. On the 
right: the corresponding residues of the momentum. 
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The two previous subsections justify the choice of a finite volume well- 
balanced scheme with a Suliciu solver adapted to treat vacuum. This scheme 
is accurate enough to study the stability of the stationary solutions and the 
asymptotic behavior of the system and we use it in the following subsections. 
In what follows, the results are displayed at time T= 100 with a space step 
equal to 0.01 and a time step satisfying the previous stability condition. 



4.3. Stability of a lateral bump. The first issue we study is the stability 
of the stationary solutions computed exactly in subsection |2j We consider 
here system ( 1.1 ) with a quadratic pressure function P{p) — sp^ and with the 
following parameters: e — D — a — h — L — 1 and x = 50. This choice guaran- 

tees that t— — — — — >0 and L>ti/\/t so the nonconstant steady state 
2eD D ' ^ ^ ^ 

with one lateral bump exists. The initial mass is taken equal to Mq = IH — . 

TT 

Assuming that the initial datum is not symmetric, the equilibrium has the 
form of the lateral bump (2.8) with the interface point 0.3943. This 



value is obtained by solving numerically (2.8c). 

We analyze the stability of this steady state under two different types 
of perturbations. First, we perturb the location of the interface point x. 



More precisely, we take initial data of the form (2.8), in which x is replaced 
by + We also recalculate the parameter the density p and the 

concentration (j) in order to have a perturbation with a zero mass. The 



parameter 5 cannot be too large, since the definition (2.8) would lead to 



negative values of the density. This type of perturbation modifies the solu- 
tion on its whole support. In the second test, we change only the density 
profile in a small region [xi,X2] satisfying < xi < X2 < x such that the initial 
density is defined as follows : 



p{x), 
p{xi), 

0, 



for < X < xi, 
for xi < X < X*, 
for X* < X < X2, 
for X2 < X < X, 
for X < X < L, 



(4.29) 



where p is the exact solution given by (2.8). The location of the jump 
X* G (xi,X2) is chosen such that the mass of the perturbation is zero. 

For the first type of perturbation with 5 = 0.1, the results are presented 
at Figure m On the left, we can see the profiles of density (on top), concen- 
tration (in the middle) and momentum (on bottom) for the initial perturbed 
data in dotted blue lines (.), the asymptotic solution after perturbation in 
solid green lines (-) and the exact solutions in lines marked with red crosses 
(+). On the right, the corresponding residues are displayed. Initially, the 
residues of density have values of order 10~^, which suggests that the solu- 
tion is still evolving. Then, at some point, they decrease rapidly to nearly 
10~^^ and stabilize at this value, which means that the solution has reached 
asymptotically the expected steady state. We remark that the asymptotic 
profiles match perfectly the expected solutions, which confirms that the sta- 



tionary solution with a lateral bump given by (2.8) is stable under this kind 
of perturbation. 
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Concentration profi 





Figure 5. On the left: Profiles for system ( [lT] ) with 7 = 2 and 
£ = D = a = b = L = l^ X = 50of density (on top) , concentration (in 
the middle) and momentum (on bottom) : (+++ red) exact so- 



lutions (2.8); (... blue) initial data with a perturbation of the in- 
terface point x*=x + 0.1; ( — green) asymptotic solutions corre- 
sponding to the previous initial data. On the right: corresponding 
residues. 



In Figu re [6| we display the results obtained with the second type of 



perturbation (4.29) with [xi,X2] = [0.6x,0.8x], with the same curves as at 
Figure [5j We notice that the numerical results confirm the stability of the 
stationary solution composed of a lateral bump. The type of perturbation 
is different, but the mechanisms of convergence to equilibrium have similar 
features. 

4.4. Dependence on the parameters x L of the asymptotic 

states. Now, in this section, we study which stationary solutions to system 



( 1.1 ) are reached asymptotically as a function of some parameters of the sys- 
tem and, in particular, how many bumps the asymptotic profile contains. 
We showed analytically that, if r < 0, the equilibrium is given by a constant 
state. The same situation occurs if t>0 and L<7v/y/r. If we increase the 
size L of the domain, nonconstant stationary solutions with several bumps 
may appear. However, we are unable to determine analytically how many 
regions of positive density the asymptotic solution will contain and a nu- 
merical study is necessary here. More precisely, our aim is to analyze how 
the chemotactic sensitivity x ^tnd the length of the domain L infiuence the 
number of bumps in the case of a quadratic pressure 7 = 2. 
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Figure 6. On the left: Profiles for system ( [lT] ) with 7 = 2 and 
£ = D = a = b = L = l^ X = 50 of density (on top) , concentration 
(in the middle) and momentum (on bottom) : (+++ red) ex- 
act solutions ( |2.8[ ); (... blue) initial data perturbed in the domain 
[^1,^2] = [0.6x,0.8x];( — green) asymptotic solutions corresponding 
to the previous initial data. On the right: corresponding residues. 

In Figure [7| we display asymptotic profiles of the density obtained for 
different lengths L of the domain. The initial data are taken equal to : 

po{x) = C (1 + sin(4^|x - L/4|)) , 00 = 0, uq^ 0, 

where ^ is computed as a function of the length L of the domain in order 
to keep the initial mass constant and equal to M = 1.3183. We observe 
that if the length of the domain is not large enough, that is L<7v/^/r^ 
the asymptotic state is constant, as expected. Increasing the size of the 
domain from L = 1 to L = 4 leads to the concentration of particles at the 
boundary and nonconstant solutions composed of one bump appear. A 
further increase, from L = 7 to L = 30, allows new bumps to appear and a 
two bumps solution is observed. 

We observe the same phenomenon when we increase the value of the 
chemosensitivity constant x- We present at Figure [8] different asymptotic 
profiles obtained for different values of x- Passing from x — ^ x — ^ leads 
to a higher concentration of particles which, as a result, produces free space 
for new bumps. However, from x = 5 to x = 200, the global form of the 
solution, i.e. two lateral bumps on each side of the domain remains the 
same, even if the two bumps tend to become higher and narrower. 
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L=5, z=3 



L = l 



L = 5 





L = 7 



L = 30 



Figure 7. Density profiles of tfie asymptotic states of system 
(1.1) with P{p)=ep'^, X = 3, e = D = a = b=l and total mass M = 
1.3183 for different values of the length of the domain L = 
{1,5,7,30}. 



X-- 



X-- 



X = 50 



X = 200 



Figure 8. Density profiles of the asymptotic states of system 



with P{p) = ep^, L = 7, £ = D = a = b = l and total mass M - 
1.3183 for different values of the chemosensitivity constant X' 
{3,5,50,200}. 



We see that the free space available for the cells has an essential effect 
on the formation of nonconstant steady states. If the size of the domain is too 
small, only constant solutions are expected. A growth of the available space, 
by increasing directly the length of the domain or by concentrating particles 
in smaller regions, leads to the formation of new bumps. However, at some 
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point, a further increase of the chemotactic force, even up to important 
values, does not produce any new bumps. So, we think that, on each domain, 
there is a maximum number of bumps once ah the parameters are fixed. 

4.5. Dependence on the adiabatic coefficient 7 of the asymptotic 
states. The adiabatic exponent 7 describes the response of cells to com- 
pression. For a high value of 7, the internal pressure repealing the cells is 
very strong. For example, for ideal gases, the value of 7 is taken between 1 
and 2, whereas in the Saint- Venant system, describing geophysical flows, it 
is equal to 2. In the case of cells, it should be much larger as they are less 
compressible than gas or water molecules. 

In Section [2j we presented a general form of nonconstant steady states 
with several bumps for arbitrary 7. However, we obtained explicit solutions 
only in the case 7 = 2. In this subsection, we study how the equilibria change 
for different values of 7 and especially how the number of bumps varies. 

We consider system ( [TT] ) with 5 = l,i:> = 0.1,a = 20,6= 10,x= 10 de- 
flned on the interval [0,3] and with initial data 

po{x) = (1.5 + sin(47r|x-L/4|)), 00 = 0, i^o = 0. 

In Figure [9| we plot the asymptotic solutions for density (on the left) and 
chemoattractant concentration (on the right) for the following values of the 
adiabatic coefficient : 7 ={2 (blue), 3 (green), 4 (red), 5 (cyan)}. We observe 
how the number of bumps changes, namely 4 bumps for 7 = 2, 3 for 7 = 3, 
2 for 7 = 4 (among which 2 lateral bumps in these 3 cases) and a constant 
profile for 7 = 5. Large values of 7 imply strong repealing forces at higher 
densities. It prevents the formation of high concentrations of cells. When 
the value of 7 increases, the height of bumps decreases, while their support 
enlarges. Moreover, when the distance between the supports of two neigh- 
boring bumps becomes zero, they join together and, for 7 high enough, the 
pressure forces are stronger than the chemotactic movement and constant 
steady states are observed. 



Density profiles at time T=200 Concentration profiles at T=200 




Figure 9. Profiles of the density p (on the left) and of the 
chemoattractant concentration (j) (on the right) at asymptotic 
states of system with £ = 1,1) = 0.1, a = 20, 6 = 10, x = 10 on the 
interval [0,3] for different values of the adiabatic exponent 7 ={2 
(blue), 3 (green), 4 (red), 5 (cyan)}. 

4.6. Dependence on the initial mass of the asymptotic states. Ex- 
periments with endothelial cells performed by Serini et.al |11| showed that 
a vascular-like network develops only if the initial density of cells ranges 
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from 100 to 400 cells/mm^. Below this interval, a disconnected structure 
is observed, while, above, a continuous carpet of cells with holes appears. 



System (1.1) was constructed to model the phenomenon of early formation 
of blood vessels and its solutions should also reflect this dependence on the 
initial mass. These experimental results suggest a particular behavior of 



steady states with several bumps of system (1.1). The regions where the 
density is strictly positive may correspond to the location of capillaries and 
would become thicker for large initial masses. Moreover, we expect to find 
a threshold value of the initial mass, above which constant equilibria are 
observed. 

We performed some simulations of system ( |1.1| ) with 7 = 2 and 7 = 3, 



which are presented at Figure llOj On the left (resp. on the right), we can 
see the asymptotic profiles for density in the case 7 = 2 (resp. 7 = 3) for 
different values of the initial mass. The initial data are the same equal to 

Po{x) = C (1 + sin(47r|x - L/4|)) , 0o = 0, uq^ 0, 

and the initial density is only multiplied by different constants in order to 
change the value of the initial mass. In the first case 7 = 2, we notice that the 
number of bumps is equal to 4 and remains exactly the same for the different 
masses. Moreover, the support of the bumps is also independent of the initial 
mass, whereas the height of the bumps increases with the initial mass. This 



can be seen theoretically, since the equations (2.8c) or (2.12c) determining 



the interface point x do not depend on the mass, while the other equations 
to determine the stationary density do. We notice that this behavior is not 
the one expected if we consider the experimental observations. However, in 
the second case 7 = 3, the dependency on the mass is totally different and 
fits the experiments mentioned above. Indeed, when the mass increases, the 
supports of the bumps become larger and the bumps join together, until 
reaching the constant equilibrium for a mass large enough. 



A 




Figure 10. Density profiles in the case P[p)=p^ (on the left) 
and P{p)=p^ (on the right) for system (1.1) with £ = 1, 1^ = 0.1, 
a = 20, 6 = 10, x = 10 and an initial datum of the form po{x) = 
^(1 + sin(47r|x — L/4|)). Comparison for different initial masses 
with ^ ={0.1 (blue), 1 (green), 5 (red), 10 (cyan)}. 
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